library(rstan)
library(rethinking)
library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
if(file.exists("~/Google Drive/Rdata/carlos.Rdata")) load("~/Google Drive/Rdata/carlos.Rdata")
tomato <- read.csv("TomatoR2CSHL.csv")
summary(tomato)
 shelf        flat            col           row            acs      trt          days           date    
 U:161   Min.   : 1.00   G      :133   Min.   :1.00   LA1954 : 40   H:495   Min.   :28.00   5/5/08:716  
 V:174   1st Qu.: 9.00   H      :127   1st Qu.:2.00   LA2695 : 39   L:513   1st Qu.:28.00   5/6/08:292  
 W:178   Median :17.00   F      :125   Median :3.00   LA1361 : 37           Median :28.00               
 X:174   Mean   :17.89   C      :117   Mean   :2.56   LA2167 : 37           Mean   :28.29               
 Y:125   3rd Qu.:28.00   D      :117   3rd Qu.:4.00   LA2773 : 37           3rd Qu.:29.00               
 Z:196   Max.   :36.00   E      :107   Max.   :4.00   LA1474 : 36           Max.   :29.00               
                         (Other):282                  (Other):782                                       
      hyp             int1            int2             int3             int4           intleng      
 Min.   : 6.17   Min.   : 0.00   Min.   : 0.000   Min.   : 0.010   Min.   : 0.030   Min.   : 0.000  
 1st Qu.:26.81   1st Qu.: 1.74   1st Qu.: 1.060   1st Qu.: 2.975   1st Qu.: 2.163   1st Qu.: 9.637  
 Median :32.02   Median : 3.59   Median : 3.120   Median : 5.625   Median : 3.995   Median :17.255  
 Mean   :33.36   Mean   : 4.71   Mean   : 4.287   Mean   : 6.794   Mean   : 5.102   Mean   :20.340  
 3rd Qu.:38.56   3rd Qu.: 6.46   3rd Qu.: 6.320   3rd Qu.: 9.367   3rd Qu.: 7.018   3rd Qu.:28.145  
 Max.   :74.60   Max.   :39.01   Max.   :28.980   Max.   :27.760   Max.   :23.280   Max.   :92.420  
                 NA's   :1       NA's   :1        NA's   :4        NA's   :102                      
    totleng          petleng         leafleng        leafwid         leafnum           ndvi    
 Min.   : 13.59   Min.   : 1.53   Min.   : 9.74   Min.   : 8.29   Min.   :3.000   Min.   :100  
 1st Qu.: 39.25   1st Qu.:11.20   1st Qu.:27.43   1st Qu.:29.48   1st Qu.:5.000   1st Qu.:108  
 Median : 50.98   Median :15.13   Median :34.59   Median :39.62   Median :5.000   Median :115  
 Mean   : 53.70   Mean   :15.92   Mean   :35.54   Mean   :39.29   Mean   :5.063   Mean   :118  
 3rd Qu.: 64.76   3rd Qu.:20.48   3rd Qu.:42.98   3rd Qu.:47.75   3rd Qu.:6.000   3rd Qu.:128  
 Max.   :129.43   Max.   :44.44   Max.   :95.19   Max.   :90.27   Max.   :8.000   Max.   :137  
                  NA's   :2       NA's   :1       NA's   :1       NA's   :1                    
      lat               lon              alt                  species      who     
 Min.   :-25.400   Min.   :-78.52   Min.   :   0   S. chilense    :207   Dan :402  
 1st Qu.:-16.607   1st Qu.:-75.92   1st Qu.:1020   S. chmielewskii:226   Pepe:606  
 Median :-14.152   Median :-73.63   Median :2240   S. habrochaites:226             
 Mean   :-14.490   Mean   :-73.71   Mean   :2035   S. pennellii   :132             
 3rd Qu.:-12.450   3rd Qu.:-71.70   3rd Qu.:3110   S. peruvianum  :217             
 Max.   : -5.767   Max.   :-68.07   Max.   :3540                                   
                                                                                   
hyp1 <- brm(hyp ~ acs*trt + (trt|species) + (1|shelf),
            prior = c(
              set_prior("normal(20,10)",class="Intercept"),
              set_prior("normal(0,10)",class="b"),
              set_prior("cauchy(0,1)",class="sigma")
            ),
            iter=5000,
            warmup=1500,
            data=tomato)
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp2 <- brm(hyp ~ acs*trt + (1|species) + (1|shelf),
            prior = c(
              set_prior("normal(20,10)",class="Intercept"),
              set_prior("normal(0,10)",class="b"),
              set_prior("cauchy(0,1)",class="sigma")
            ),
            iter=5000,
            warmup=1500,
            data=tomato)
Compiling the C++ model
Start sampling
starting worker pid=99011 on localhost:11828 at 12:22:36.419
starting worker pid=99019 on localhost:11828 at 12:22:36.592
starting worker pid=99027 on localhost:11828 at 12:22:36.760
starting worker pid=99035 on localhost:11828 at 12:22:36.935

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 1501 / 5000 [ 30%]  (Sampling)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1501 / 5000 [ 30%]  (Sampling)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1501 / 5000 [ 30%]  (Sampling)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 1501 / 5000 [ 30%]  (Sampling)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 180.896 seconds (Warm-up)
               226.983 seconds (Sampling)
               407.879 seconds (Total)


Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 180.073 seconds (Warm-up)
               230.852 seconds (Sampling)
               410.924 seconds (Total)


Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 208.94 seconds (Warm-up)
               243.667 seconds (Sampling)
               452.607 seconds (Total)


Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 171.104 seconds (Warm-up)
               349.297 seconds (Sampling)
               520.401 seconds (Total)
There were 541 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp3 <- update(hyp1, formula. = hyp ~ acs*trt + (1|species))
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=99334 on localhost:11828 at 12:32:57.610
starting worker pid=99342 on localhost:11828 at 12:32:57.791
starting worker pid=99350 on localhost:11828 at 12:32:57.969
starting worker pid=99358 on localhost:11828 at 12:32:58.153

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 195.425 seconds (Warm-up)
               63.3163 seconds (Sampling)
               258.741 seconds (Total)


Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 199.978 seconds (Warm-up)
               107.182 seconds (Sampling)
               307.16 seconds (Total)


Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 205.032 seconds (Warm-up)
               114.488 seconds (Sampling)
               319.52 seconds (Total)


Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 207.932 seconds (Warm-up)
               119.223 seconds (Sampling)
               327.155 seconds (Total)
There were 1105 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp4 <- update(hyp1, formula. = hyp ~ acs*trt + (1|shelf))
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=99586 on localhost:11828 at 12:39:50.457
starting worker pid=99594 on localhost:11828 at 12:39:50.613
starting worker pid=99602 on localhost:11828 at 12:39:50.770
starting worker pid=99610 on localhost:11828 at 12:39:50.922

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 151.259 seconds (Warm-up)
               78.1472 seconds (Sampling)
               229.407 seconds (Total)


Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 138.363 seconds (Warm-up)
               105.666 seconds (Sampling)
               244.029 seconds (Total)


Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 153.591 seconds (Warm-up)
               94.0214 seconds (Sampling)
               247.613 seconds (Total)


Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 144.983 seconds (Warm-up)
               136.555 seconds (Sampling)
               281.538 seconds (Total)
There were 7 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp5 <- update(hyp1, formula. = hyp ~ acs*trt)
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=99715 on localhost:11828 at 12:45:40.313
starting worker pid=99723 on localhost:11828 at 12:45:40.486
starting worker pid=99731 on localhost:11828 at 12:45:40.654
starting worker pid=99739 on localhost:11828 at 12:45:40.814

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 85.807 seconds (Warm-up)
               39.2256 seconds (Sampling)
               125.033 seconds (Total)


Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 86.8877 seconds (Warm-up)
               36.1756 seconds (Sampling)
               123.063 seconds (Total)


Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 91.221 seconds (Warm-up)
               36.5213 seconds (Sampling)
               127.742 seconds (Total)


Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 90.6127 seconds (Warm-up)
               38.1106 seconds (Sampling)
               128.723 seconds (Total)
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
loo(hyp1,hyp2,hyp3,hyp4,hyp5)
Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
              LOOIC    SE
hyp1        7042.25 58.38
hyp2        7042.55 58.24
hyp3        7088.25 58.10
hyp4        7041.43 58.26
hyp5        7084.61 58.15
hyp1 - hyp2   -0.30  1.00
hyp1 - hyp3  -46.00 12.92
hyp1 - hyp4    0.82  1.07
hyp1 - hyp5  -42.36 12.85
hyp2 - hyp3  -45.70 12.78
hyp2 - hyp4    1.11  0.48
hyp2 - hyp5  -42.06 12.71
hyp3 - hyp4   46.81 12.76
hyp3 - hyp5    3.64  0.84
hyp4 - hyp5  -43.18 12.68

So we should go with hyp4 (no species effect)

summary(hyp4)
 Family: gaussian (identity) 
Formula: hyp ~ acs + trt + (1 | shelf) + acs:trt 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
   WAIC: Not computed
 
Group-Level Effects: 
~shelf (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     2.85      1.34     1.23     6.36       3387    1

Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept         28.83      2.09    24.57    32.97       3863    1
acsLA1305         -0.05      2.45    -4.96     4.65       4684    1
acsLA1306         -1.56      2.22    -5.91     2.82       5142    1
acsLA1316          9.38      2.24     5.01    13.74       4346    1
acsLA1317          5.73      2.16     1.51     9.92       4888    1
acsLA1318         -1.88      2.19    -6.23     2.47       4024    1
acsLA1336         11.36      2.29     6.81    15.83       4572    1
acsLA1361          1.46      2.07    -2.76     5.50       4150    1
acsLA1363         -0.67      2.16    -4.92     3.57       4379    1
acsLA1367         -8.24      2.81   -13.80    -2.75      10000    1
acsLA1474         12.64      2.13     8.49    16.80       4482    1
acsLA1691         -1.79      2.31    -6.34     2.76       4669    1
acsLA1912         -5.45      2.57   -10.55    -0.45       5795    1
acsLA1918          6.63      2.30     2.15    11.22       5424    1
acsLA1928          0.71      3.97    -6.95     8.65      10000    1
acsLA1954          7.32      2.05     3.34    11.33       3930    1
acsLA1958          7.44      2.16     3.24    11.69       4442    1
acsLA1969         -4.93      2.26    -9.37    -0.49       4521    1
acsLA1973          6.00      2.32     1.48    10.42       5604    1
acsLA2167         -1.29      2.08    -5.38     2.78       3931    1
acsLA2204         -1.73      2.05    -5.76     2.24       3863    1
acsLA2560         -5.71      3.88   -13.30     1.91      10000    1
acsLA2580         -6.53      2.67   -11.78    -1.28       5300    1
acsLA2663         -2.22      2.14    -6.40     1.96       4201    1
acsLA2695          0.31      2.08    -3.76     4.37       4053    1
acsLA2748         -1.61      2.56    -6.60     3.37       5837    1
acsLA2773          4.64      2.09     0.59     8.79       4160    1
acsLA2884          1.25      2.34    -3.31     5.81       4433    1
acsLA2891          0.65      2.38    -4.08     5.24       5447    1
acsLA2931          8.31      2.30     3.81    12.82       4775    1
acsLA361           2.72      2.36    -1.85     7.34       4959    1
acsLA3788         -4.77      2.27    -9.19    -0.32       4992    1
acsLA3791         -1.46      2.45    -6.24     3.32       5022    1
acsLA3795          0.67      2.24    -3.73     4.98       4991    1
acsLA446           1.57      2.36    -3.08     6.17       4771    1
acsLA716          -9.74      3.46   -16.41    -2.87      10000    1
trtL               3.66      2.76    -2.09     9.09       3947    1
acsLA1305:trtL    -5.66      3.16   -11.85     0.45       6853    1
acsLA1306:trtL    -1.23      2.86    -6.94     4.36       6470    1
acsLA1316:trtL     0.88      2.90    -4.76     6.48       6674    1
acsLA1317:trtL    -1.38      2.82    -6.92     4.17      10000    1
acsLA1318:trtL    -0.44      2.95    -6.16     5.45       6977    1
acsLA1336:trtL     3.79      2.93    -1.95     9.53       6621    1
acsLA1361:trtL     2.74      2.78    -2.77     8.22      10000    1
acsLA1363:trtL    -2.44      2.86    -7.97     3.25       6157    1
acsLA1367:trtL     9.39      3.78     2.00    16.80      10000    1
acsLA1474:trtL     3.83      2.80    -1.69     9.15      10000    1
acsLA1691:trtL    -0.61      3.14    -6.77     5.48      10000    1
acsLA1912:trtL     0.30      3.26    -6.07     6.52       6838    1
acsLA1918:trtL    -1.63      3.26    -8.04     4.82      10000    1
acsLA1928:trtL    -9.85      5.75   -21.04     1.63      10000    1
acsLA1954:trtL     1.68      2.68    -3.58     6.96       6278    1
acsLA1958:trtL     2.17      2.89    -3.56     7.81      10000    1
acsLA1969:trtL    -0.88      3.08    -6.85     5.11       5883    1
acsLA1973:trtL     8.11      3.02     2.29    14.08      10000    1
acsLA2167:trtL    -0.54      2.78    -6.07     4.96       5263    1
acsLA2204:trtL     1.20      2.83    -4.48     6.68      10000    1
acsLA2560:trtL     6.76      4.80    -2.62    16.17      10000    1
acsLA2580:trtL     4.17      3.38    -2.30    10.85      10000    1
acsLA2663:trtL    -1.37      2.81    -6.96     4.08       5960    1
acsLA2695:trtL    -0.78      2.74    -6.11     4.60       6219    1
acsLA2748:trtL    -0.40      3.30    -6.95     6.13      10000    1
acsLA2773:trtL     0.09      2.80    -5.54     5.44       6499    1
acsLA2884:trtL     1.70      3.07    -4.40     7.64      10000    1
acsLA2891:trtL     7.06      3.20     0.91    13.22      10000    1
acsLA2931:trtL     4.77      2.94    -0.98    10.42      10000    1
acsLA361:trtL      0.80      3.05    -5.03     6.76      10000    1
acsLA3788:trtL    10.09      3.21     3.73    16.45      10000    1
acsLA3791:trtL     3.37      3.09    -2.80     9.39      10000    1
acsLA3795:trtL     3.31      2.95    -2.51     9.12      10000    1
acsLA446:trtL     -0.47      3.27    -6.88     5.88      10000    1
acsLA716:trtL      4.76      4.35    -3.81    13.07      10000    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     7.67      0.18     7.34     8.03      10000    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(hyp4, ask=FALSE)

ranef(hyp4)
$shelf
   Intercept
U -1.4409542
V  1.7748644
W  0.6353095
X -2.0373836
Y -0.1921667
Z  2.4487958
pred.df <- unique(tomato[,c("species","acs","trt")])
head(pred.df)
hyp.prediction <- cbind(pred.df,fitted(hyp4,newdata = pred.df, re_formula = NA))
colnames(hyp.prediction)[6:7] <- c("low95","high95")
head(hyp.prediction)
write.csv(hyp.prediction,file="hyp_accession.csv")
pl <- ggplot(hyp.prediction,aes(x=acs,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl <- pl + facet_wrap( ~ species, ncol=2, scales="free_x")
pl + ggtitle("hypocotyls per accession")
ggsave("hypocotyls_accession.pdf",height=11,width=8.5)

by(tomato$intleng,tomato$trt,mean)
int1 <- brm(intleng ~ acs*trt + (trt|species) + (1|shelf),
            prior = c(
              set_prior("normal(10,10)",class="Intercept"),
              set_prior("normal(0,20)",class="b"),
              set_prior("cauchy(0,1)",class="sigma")
            ),
            iter=5000,
            warmup=1500,
            data=tomato)
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
int2 <- update(int1, formula. = intleng ~ acs*trt + (1|species) + (1|shelf))
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=177 on localhost:11828 at 13:03:31.512
starting worker pid=185 on localhost:11828 at 13:03:31.684
starting worker pid=197 on localhost:11828 at 13:03:31.854
starting worker pid=205 on localhost:11828 at 13:03:32.026

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 228.748 seconds (Warm-up)
               165.635 seconds (Sampling)
               394.382 seconds (Total)


Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 238.932 seconds (Warm-up)
               163.406 seconds (Sampling)
               402.338 seconds (Total)


Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 258.535 seconds (Warm-up)
               149.635 seconds (Sampling)
               408.17 seconds (Total)


Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 229.139 seconds (Warm-up)
               265.832 seconds (Sampling)
               494.97 seconds (Total)
There were 533 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
int3 <- update(int2, formula. = intleng ~ acs*trt + (1|species))
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=502 on localhost:11828 at 13:16:43.400
starting worker pid=515 on localhost:11828 at 13:16:43.568
starting worker pid=524 on localhost:11828 at 13:16:43.724
starting worker pid=532 on localhost:11828 at 13:16:43.880

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 222.004 seconds (Warm-up)
               31.0146 seconds (Sampling)
               253.019 seconds (Total)


Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 206.206 seconds (Warm-up)
               76.8814 seconds (Sampling)
               283.087 seconds (Total)


Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 212.625 seconds (Warm-up)
               117.46 seconds (Sampling)
               330.086 seconds (Total)


Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 240.151 seconds (Warm-up)
               95.9457 seconds (Sampling)
               336.097 seconds (Total)
There were 2407 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
int4 <- update(int3, formula. = intleng ~ acs*trt + (1|shelf))
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=838 on localhost:11828 at 13:27:07.734
starting worker pid=848 on localhost:11828 at 13:27:07.915
starting worker pid=863 on localhost:11828 at 13:27:08.097
starting worker pid=875 on localhost:11828 at 13:27:08.272

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 126.289 seconds (Warm-up)
               70.5382 seconds (Sampling)
               196.827 seconds (Total)


Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 133.997 seconds (Warm-up)
               71.3704 seconds (Sampling)
               205.367 seconds (Total)


Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 129.769 seconds (Warm-up)
               79.0085 seconds (Sampling)
               208.778 seconds (Total)


Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 138.738 seconds (Warm-up)
               78.9029 seconds (Sampling)
               217.641 seconds (Total)
There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupExamine the pairs() plot to diagnose sampling problems
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
int5 <- update(int1, formula. = intleng ~ acs*trt)
The desired formula changes require recompling the model
Compiling the C++ model
Start sampling
starting worker pid=1172 on localhost:11828 at 13:32:53.126
starting worker pid=1180 on localhost:11828 at 13:32:53.309
starting worker pid=1188 on localhost:11828 at 13:32:53.487
starting worker pid=1196 on localhost:11828 at 13:32:53.668

SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 1).

Chain 1, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 2).

Chain 2, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 3).

Chain 3, Iteration:    1 / 5000 [  0%]  (Warmup)
SAMPLING FOR MODEL 'gaussian(identity) brms-model' NOW (CHAIN 4).

Chain 4, Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4, Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4, Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4, Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 2, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 4, Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4, Iteration: 2501 / 5000 [ 50%]  (Sampling)
Chain 3, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 4, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2, Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2, Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 90.9246 seconds (Warm-up)
               44.8574 seconds (Sampling)
               135.782 seconds (Total)


Chain 2, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1, Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 4, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 97.5076 seconds (Warm-up)
               49.4388 seconds (Sampling)
               146.946 seconds (Total)


Chain 2, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1, Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 96.0257 seconds (Warm-up)
               62.106 seconds (Sampling)
               158.132 seconds (Total)


Chain 1, Iteration: 5000 / 5000 [100%]  (Sampling)
 Elapsed Time: 86.7937 seconds (Warm-up)
               75.9581 seconds (Sampling)
               162.752 seconds (Total)
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
Warning messages:
1: package ‘rstan’ was built under R version 3.3.2 
2: package ‘ggplot2’ was built under R version 3.3.2 
3: package ‘StanHeaders’ was built under R version 3.3.2 
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
loo(int1,int2,int3,int4,int5)
Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
              LOOIC    SE
int1        7396.58 64.57
int2        7399.29 64.61
int3        7413.85 62.36
int4        7397.58 64.51
int5        7412.23 62.83
int1 - int2   -2.71  0.92
int1 - int3  -17.28 10.76
int1 - int4   -1.00  0.82
int1 - int5  -15.66 10.50
int2 - int3  -14.56 10.70
int2 - int4    1.71  0.65
int2 - int5  -12.94 10.44
int3 - int4   16.28 10.59
int3 - int5    1.62  2.20
int4 - int5  -14.66 10.33

So we should go with int4 (no species effect). Case could also be made for int5 (no random effects)

summary(int4)
 Family: gaussian (identity) 
Formula: intleng ~ acs + trt + (1 | shelf) + acs:trt 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
   WAIC: Not computed
 
Group-Level Effects: 
~shelf (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     2.35      1.33     0.82     5.72       2915    1

Population-Level Effects: 
               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept          8.80      2.48     3.79    13.60       1058    1
acsLA1305          3.38      3.29    -2.98     9.89       1902    1
acsLA1306          1.99      3.02    -3.92     7.70       1618    1
acsLA1316          0.79      3.04    -5.17     6.75       1603    1
acsLA1317         11.85      3.00     5.97    17.75       1708    1
acsLA1318          7.79      3.00     1.83    13.57       1649    1
acsLA1336          7.23      3.10     1.08    13.32       1678    1
acsLA1361          5.39      2.83    -0.17    10.96       1542    1
acsLA1363         11.02      2.93     5.14    16.80       1593    1
acsLA1367          0.15      3.65    -6.99     7.42       2263    1
acsLA1474          1.50      2.87    -4.21     7.08       1518    1
acsLA1691          1.20      3.10    -4.85     7.20       1665    1
acsLA1912         -2.30      3.47    -9.10     4.49       2375    1
acsLA1918          3.29      3.14    -2.92     9.55       1739    1
acsLA1928          2.26      5.24    -7.92    12.45       4792    1
acsLA1954          6.35      2.82     0.83    11.87       1457    1
acsLA1958         -1.35      2.96    -7.15     4.47       1670    1
acsLA1969         -6.59      3.06   -12.63    -0.54       1628    1
acsLA1973         11.55      3.10     5.50    17.64       1627    1
acsLA2167          3.58      2.87    -2.08     9.14       1573    1
acsLA2204          8.46      2.83     2.92    13.96       1368    1
acsLA2560          0.13      5.20   -10.02    10.33       4298    1
acsLA2580         -4.68      3.56   -11.68     2.28       2107    1
acsLA2663          4.95      2.96    -0.88    10.74       1573    1
acsLA2695          3.41      2.87    -2.33     8.95       1427    1
acsLA2748         -2.50      3.48    -9.50     4.35       2026    1
acsLA2773         -3.84      2.87    -9.61     1.74       1532    1
acsLA2884         -2.37      3.08    -8.39     3.65       1658    1
acsLA2891         -2.03      3.18    -8.30     4.17       1804    1
acsLA2931         -3.28      3.08    -9.27     2.79       1627    1
acsLA361          12.37      3.20     5.92    18.58       1791    1
acsLA3788         -3.55      3.10    -9.65     2.65       1627    1
acsLA3791         -6.13      3.27   -12.57     0.29       1837    1
acsLA3795         -0.67      3.04    -6.65     5.26       1604    1
acsLA446           2.15      3.11    -4.02     8.18       1657    1
acsLA716          -2.32      4.71   -11.61     6.92       3806    1
trtL              12.68      3.25     6.06    19.04       1177    1
acsLA1305:trtL     1.89      4.22    -6.43    10.14       2225    1
acsLA1306:trtL     5.07      3.90    -2.55    12.83       1968    1
acsLA1316:trtL    -1.84      3.96    -9.49     6.02       1867    1
acsLA1317:trtL     7.93      3.95     0.40    15.77       1736    1
acsLA1318:trtL     6.87      4.05    -1.05    14.91       2048    1
acsLA1336:trtL     7.23      4.07    -0.80    15.20       1797    1
acsLA1361:trtL    15.69      3.85     8.21    23.35       1884    1
acsLA1363:trtL     5.65      3.90    -1.89    13.30       1866    1
acsLA1367:trtL     8.91      4.97    -0.95    18.60       3153    1
acsLA1474:trtL     0.51      3.82    -6.87     8.09       1778    1
acsLA1691:trtL    -5.80      4.22   -14.15     2.39       2193    1
acsLA1912:trtL     2.23      4.44    -6.38    10.99       2420    1
acsLA1918:trtL    -5.97      4.44   -14.55     2.65       2279    1
acsLA1928:trtL    -1.75      7.81   -17.15    13.58      10000    1
acsLA1954:trtL     1.84      3.72    -5.36     9.26       1746    1
acsLA1958:trtL     6.58      3.96    -1.24    14.34       2045    1
acsLA1969:trtL     6.55      4.15    -1.56    14.84       2226    1
acsLA1973:trtL    15.81      4.05     7.95    23.79       2075    1
acsLA2167:trtL     0.88      3.84    -6.64     8.41       1815    1
acsLA2204:trtL     5.15      3.90    -2.41    12.78       1836    1
acsLA2560:trtL    10.77      6.49    -1.93    23.36       4647    1
acsLA2580:trtL     8.10      4.60    -0.83    17.04       2313    1
acsLA2663:trtL    -0.65      3.90    -8.23     7.11       1835    1
acsLA2695:trtL    -4.27      3.76   -11.62     3.19       1683    1
acsLA2748:trtL    17.95      4.54     8.96    26.81       2430    1
acsLA2773:trtL    13.28      3.83     5.76    20.82       1765    1
acsLA2884:trtL     0.50      4.09    -7.45     8.55       1936    1
acsLA2891:trtL     8.71      4.23     0.34    16.92       2154    1
acsLA2931:trtL     9.45      3.96     1.67    17.27       1823    1
acsLA361:trtL      8.33      4.20     0.15    16.69       2112    1
acsLA3788:trtL     4.89      4.29    -3.43    13.33       2042    1
acsLA3791:trtL     1.49      4.15    -6.51     9.71       2233    1
acsLA3795:trtL     4.81      4.02    -3.23    12.67       1841    1
acsLA446:trtL      6.34      4.32    -2.07    14.82       2257    1
acsLA716:trtL     -0.28      6.05   -12.14    11.42       3960    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     9.12      0.21     8.72     9.54      10000    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(int4, ask=FALSE)

ranef(int4)
$shelf
    Intercept
U -1.60562115
V  2.34342216
W -0.27581150
X -0.01534482
Y -0.56118509
Z  0.73463441
intleng.prediction <- cbind(pred.df,fitted(int4,newdata = pred.df, re_formula = NA))
colnames(intleng.prediction)[6:7] <- c("low95","high95")
head(intleng.prediction)
write.csv(intleng.prediction,file="intleng_accession.csv")
pl <- ggplot(intleng.prediction,aes(x=acs,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl <- pl + facet_wrap( ~ species, ncol=2, scales="free_x")
pl + ggtitle("internode")
ggsave("internodes_accession.pdf",height=11,width=8.5)

species level analyses

hypocotyl

hyp6 <- brm(hyp ~ species*trt + (trt|acs) + (1|shelf),
            prior = c(
              set_prior("normal(20,10)",class="Intercept"),
              set_prior("normal(0,10)",class="b"),
              set_prior("cauchy(0,1)",class="sigma")
            ),
            iter=5000,
            warmup=1500,
            data=tomato)
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp7 <- update(hyp6, formula. = . ~ species*trt + (1|acs) + (1|shelf) )
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp8 <- update(hyp7, formula. = . ~ species*trt + (1|shelf) )
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp9 <- update(hyp7, formula. = . ~ species*trt + (1|acs) )
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
hyp10 <- update(hyp7, formula. = . ~ species*trt)
save.image(file="~/Google Drive/Rdata/carlos.Rdata")
summary(hyp6)
 Family: gaussian (identity) 
Formula: hyp ~ species * trt + (trt | acs) + (1 | shelf) 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 1500; thin = 1; 
         total post-warmup samples = 14000
   WAIC: Not computed
 
Group-Level Effects: 
~acs (Number of levels: 36) 
                    Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)           4.15      0.69     2.99     5.68       5208    1
sd(trtL)                1.81      0.71     0.45     3.28       4807    1
cor(Intercept,trtL)     0.67      0.27     0.01     0.99       6552    1

~shelf (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     2.88      1.37     1.24     6.49       4139    1

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                     30.90      2.46    25.93    35.65       3643    1
speciesS.chmielewskii         -0.51      2.38    -5.14     4.19       4060    1
speciesS.habrochaites         -1.60      2.32    -6.17     3.01       3681    1
speciesS.pennellii            -8.41      2.49   -13.26    -3.47       4310    1
speciesS.peruvianum            3.52      2.40    -1.19     8.24       3671    1
trtL                           5.50      2.68    -0.07    10.59       4802    1
speciesS.chmielewskii:trtL    -3.04      1.72    -6.44     0.33       8131    1
speciesS.habrochaites:trtL    -2.30      1.71    -5.72     1.09       8232    1
speciesS.pennellii:trtL        4.31      1.93     0.46     8.12       8650    1
speciesS.peruvianum:trtL       0.19      1.73    -3.19     3.60       7729    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     7.67      0.18     7.33     8.02      14000    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
summary(hyp7)
 Family: gaussian (identity) 
Formula: hyp ~ species + trt + (1 | acs) + (1 | shelf) + species:trt 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
   WAIC: Not computed
 
Group-Level Effects: 
~acs (Number of levels: 36) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     4.85      0.69     3.69      6.4       2597    1

~shelf (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     3.02       1.7     1.22     7.53        627 1.01

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                     30.83      2.54    25.74    35.77       3164    1
speciesS.chmielewskii         -0.48      2.63    -5.66     4.78       3295    1
speciesS.habrochaites         -1.44      2.50    -6.35     3.33       3365    1
speciesS.pennellii            -8.34      2.66   -13.58    -3.04       3809    1
speciesS.peruvianum            3.49      2.57    -1.56     8.53       3504    1
trtL                           5.56      2.74    -0.34    10.83       1666    1
speciesS.chmielewskii:trtL    -3.19      1.45    -6.02    -0.33       6929    1
speciesS.habrochaites:trtL    -2.38      1.45    -5.23     0.45       4668    1
speciesS.pennellii:trtL        4.17      1.64     0.99     7.44       6833    1
speciesS.peruvianum:trtL       0.26      1.44    -2.63     3.04       7030    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     7.71      0.18     7.37     8.06      10000    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
summary(hyp8)
 Family: gaussian (identity) 
Formula: hyp ~ species + trt + (1 | shelf) + species:trt 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
   WAIC: Not computed
 
Group-Level Effects: 
~shelf (Number of levels: 6) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     2.53      1.23     1.01     5.78       2411    1

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                     31.47      1.70    27.97    34.90       3964    1
speciesS.chmielewskii         -1.06      1.18    -3.38     1.26       6018    1
speciesS.habrochaites         -2.06      1.14    -4.29     0.13       5576    1
speciesS.pennellii            -8.47      1.38   -11.20    -5.81       6226    1
speciesS.peruvianum            3.66      1.18     1.32     5.98       5608    1
trtL                           6.01      2.40     1.11    10.73       3406    1
speciesS.chmielewskii:trtL    -3.23      1.63    -6.44    -0.04       5875    1
speciesS.habrochaites:trtL    -2.67      1.62    -5.85     0.45       5471    1
speciesS.pennellii:trtL        3.49      1.86    -0.10     7.16       6458    1
speciesS.peruvianum:trtL       0.03      1.63    -3.14     3.22       5323    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     8.88       0.2      8.5     9.27      10000    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
summary(hyp9)
 Family: gaussian (identity) 
Formula: hyp ~ species + trt + (1 | acs) + species:trt 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
   WAIC: Not computed
 
Group-Level Effects: 
~acs (Number of levels: 36) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     4.74      0.69      3.6     6.29       2491    1

Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                     30.98      1.81    27.39    34.58       1700    1
speciesS.chmielewskii         -0.44      2.57    -5.54     4.59       2024    1
speciesS.habrochaites         -1.50      2.50    -6.40     3.47       1920    1
speciesS.pennellii            -7.64      2.68   -13.02    -2.38       2039    1
speciesS.peruvianum            3.51      2.59    -1.68     8.53       1903    1
trtL                           6.07      1.05     4.02     8.13       2498    1
speciesS.chmielewskii:trtL    -3.33      1.48    -6.22    -0.47       3381    1
speciesS.habrochaites:trtL    -2.62      1.46    -5.49     0.23       3241    1
speciesS.pennellii:trtL        3.58      1.71     0.23     6.91       4482    1
speciesS.peruvianum:trtL       0.20      1.49    -2.72     3.07       3708    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma      7.9      0.18     7.56     8.26       9577    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
summary(hyp10)
 Family: gaussian (identity) 
Formula: hyp ~ species + trt + species:trt 
   Data: tomato (Number of observations: 1008) 
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1; 
         total post-warmup samples = 10000
   WAIC: Not computed
 
Population-Level Effects: 
                           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept                     31.53      0.82    29.92    33.15       4801    1
speciesS.chmielewskii         -0.98      1.15    -3.20     1.27       5388    1
speciesS.habrochaites         -2.07      1.12    -4.28     0.12       5441    1
speciesS.pennellii            -7.82      1.34   -10.42    -5.20       5804    1
speciesS.peruvianum            3.72      1.17     1.37     6.03       5668    1
trtL                           6.28      1.14     4.05     8.53       4160    1
speciesS.chmielewskii:trtL    -3.32      1.60    -6.48    -0.24       4773    1
speciesS.habrochaites:trtL    -2.86      1.59    -5.93     0.23       4837    1
speciesS.pennellii:trtL        3.08      1.83    -0.53     6.66       5175    1
speciesS.peruvianum:trtL       0.00      1.62    -3.16     3.19       5134    1

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma        9      0.19     8.63     9.39       8766    1

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
loo(hyp6,hyp7,hyp8,hyp9,hyp10)
               LOOIC    SE
hyp6         7019.05 58.82
hyp7         7022.09 58.24
hyp8         7277.38 58.08
hyp9         7068.53 58.00
hyp10        7300.36 57.60
hyp6 - hyp7    -3.04  6.03
hyp6 - hyp8  -258.33 31.40
hyp6 - hyp9   -49.48 14.03
hyp6 - hyp10 -281.31 32.13
hyp7 - hyp8  -255.29 30.26
hyp7 - hyp9   -46.45 12.72
hyp7 - hyp10 -278.27 31.04
hyp8 - hyp9   208.85 31.60
hyp8 - hyp10  -22.98  9.62
hyp9 - hyp10 -231.83 28.97

hyp6, 7 equivalent (no need for trt|acs) hyp8, 10 much worse (keep acs intercept random effect) hyp9 a bit worse (some evidenc for shelf effects)

go with hyp7

pred.df.species <- unique(tomato[,c("species","trt")])
hyp.species.prediction <- cbind(pred.df.species,fitted(hyp7,newdata = pred.df.species, re_formula = NA))
colnames(hyp.species.prediction)[5:6] <- c("low95","high95")
head(hyp.species.prediction)
write.csv(hyp.species.prediction,file="hyp_species.csv")
pl <- ggplot(hyp.species.prediction,aes(x=species,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl + ggtitle("hypocotyls per species")
ggsave("hypocotyls_species.pdf")
Saving 6.77 x 4.18 in image

internode

loo(int6,int7,int8,int9,int10)
Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
               LOOIC    SE
int6         7371.26 65.23
int7         7419.36 64.33
int8         7635.99 66.00
int9         7433.97 62.54
int10        7644.33 64.68
int6 - int7   -48.10 17.12
int6 - int8  -264.73 38.16
int6 - int9   -62.71 19.68
int6 - int10 -273.08 38.97
int7 - int8  -216.63 27.97
int7 - int9   -14.61  9.84
int7 - int10 -224.97 29.25
int8 - int9   202.02 28.99
int8 - int10   -8.34  8.45
int9 - int10 -210.37 27.57

int6 is best

int.species.prediction <- cbind(pred.df.species,fitted(int6,newdata = pred.df.species, re_formula = NA))
colnames(int.species.prediction)[5:6] <- c("low95","high95")
head(int.species.prediction)
write.csv(int.species.prediction,file="int_species.csv")
pl <- ggplot(int.species.prediction,aes(x=species,y=Estimate,ymin=low95,ymax=high95, fill=trt))
pl <- pl + geom_bar(position="dodge",stat="identity")
pl <- pl + geom_errorbar(position=position_dodge(width=.9), width=0.5)
pl + ggtitle("internode per species")
ggsave("internodes_species.pdf")
Saving 6.77 x 4.18 in image

LS0tCnRpdGxlOiAiVG9tYXRvIGZvciBDYXJsb3MiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0gCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCBjYWNoZSA9IFRSVUUsIGF1dG9kZXAgPSBUUlVFKQpgYGAKCmBgYHtyfQpsaWJyYXJ5KHJzdGFuKQpsaWJyYXJ5KHJldGhpbmtpbmcpCmxpYnJhcnkoYnJtcykKcnN0YW5fb3B0aW9ucyhhdXRvX3dyaXRlID0gVFJVRSkKb3B0aW9ucyhtYy5jb3JlcyA9IHBhcmFsbGVsOjpkZXRlY3RDb3JlcygpKQppZihmaWxlLmV4aXN0cygifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikpIGxvYWQoIn4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CnRvbWF0byA8LSByZWFkLmNzdigiVG9tYXRvUjJDU0hMLmNzdiIpCnN1bW1hcnkodG9tYXRvKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwMSA8LSBicm0oaHlwIH4gYWNzKnRydCArICh0cnR8c3BlY2llcykgKyAoMXxzaGVsZiksCiAgICAgICAgICAgIHByaW9yID0gYygKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgyMCwxMCkiLGNsYXNzPSJJbnRlcmNlcHQiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEwKSIsY2xhc3M9ImIiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoImNhdWNoeSgwLDEpIixjbGFzcz0ic2lnbWEiKQogICAgICAgICAgICApLAogICAgICAgICAgICBpdGVyPTUwMDAsCiAgICAgICAgICAgIHdhcm11cD0xNTAwLAogICAgICAgICAgICBkYXRhPXRvbWF0bykKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQpoeXAyIDwtIGJybShoeXAgfiBhY3MqdHJ0ICsgKDF8c3BlY2llcykgKyAoMXxzaGVsZiksCiAgICAgICAgICAgIHByaW9yID0gYygKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgyMCwxMCkiLGNsYXNzPSJJbnRlcmNlcHQiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEwKSIsY2xhc3M9ImIiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoImNhdWNoeSgwLDEpIixjbGFzcz0ic2lnbWEiKQogICAgICAgICAgICApLAogICAgICAgICAgICBpdGVyPTUwMDAsCiAgICAgICAgICAgIHdhcm11cD0xNTAwLAogICAgICAgICAgICBkYXRhPXRvbWF0bykKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQpoeXAzIDwtIHVwZGF0ZShoeXAxLCBmb3JtdWxhLiA9IGh5cCB+IGFjcyp0cnQgKyAoMXxzcGVjaWVzKSkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQpoeXA0IDwtIHVwZGF0ZShoeXAxLCBmb3JtdWxhLiA9IGh5cCB+IGFjcyp0cnQgKyAoMXxzaGVsZikpCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KaHlwNSA8LSB1cGRhdGUoaHlwMSwgZm9ybXVsYS4gPSBoeXAgfiBhY3MqdHJ0KQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CmxvbyhoeXAxLGh5cDIsaHlwMyxoeXA0LGh5cDUpCmBgYAoKU28gd2Ugc2hvdWxkIGdvIHdpdGggaHlwNCAobm8gc3BlY2llcyBlZmZlY3QpCmBgYHtyfQpzdW1tYXJ5KGh5cDQpCnBsb3QoaHlwNCwgYXNrPUZBTFNFKQpyYW5lZihoeXA0KQpgYGAKCmBgYHtyfQpwcmVkLmRmIDwtIHVuaXF1ZSh0b21hdG9bLGMoInNwZWNpZXMiLCJhY3MiLCJ0cnQiKV0pCmhlYWQocHJlZC5kZikKaHlwLnByZWRpY3Rpb24gPC0gY2JpbmQocHJlZC5kZixmaXR0ZWQoaHlwNCxuZXdkYXRhID0gcHJlZC5kZiwgcmVfZm9ybXVsYSA9IE5BKSkKY29sbmFtZXMoaHlwLnByZWRpY3Rpb24pWzY6N10gPC0gYygibG93OTUiLCJoaWdoOTUiKQpoZWFkKGh5cC5wcmVkaWN0aW9uKQp3cml0ZS5jc3YoaHlwLnByZWRpY3Rpb24sZmlsZT0iaHlwX2FjY2Vzc2lvbi5jc3YiKQpgYGAKCmBgYHtyfQpwbCA8LSBnZ3Bsb3QoaHlwLnByZWRpY3Rpb24sYWVzKHg9YWNzLHk9RXN0aW1hdGUseW1pbj1sb3c5NSx5bWF4PWhpZ2g5NSwgZmlsbD10cnQpKQpwbCA8LSBwbCArIGdlb21fYmFyKHBvc2l0aW9uPSJkb2RnZSIsc3RhdD0iaWRlbnRpdHkiKQpwbCA8LSBwbCArIGdlb21fZXJyb3JiYXIocG9zaXRpb249cG9zaXRpb25fZG9kZ2Uod2lkdGg9LjkpLCB3aWR0aD0wLjUpCnBsIDwtIHBsICsgZmFjZXRfd3JhcCggfiBzcGVjaWVzLCBuY29sPTIsIHNjYWxlcz0iZnJlZV94IikKcGwgKyBnZ3RpdGxlKCJoeXBvY290eWxzIHBlciBhY2Nlc3Npb24iKQpnZ3NhdmUoImh5cG9jb3R5bHNfYWNjZXNzaW9uLnBkZiIsaGVpZ2h0PTExLHdpZHRoPTguNSkKYGBgCgpgYGB7ciwgcmVzdWx0cz0naGlkZSd9CmJ5KHRvbWF0byRpbnRsZW5nLHRvbWF0byR0cnQsbWVhbikKaW50MSA8LSBicm0oaW50bGVuZyB+IGFjcyp0cnQgKyAodHJ0fHNwZWNpZXMpICsgKDF8c2hlbGYpLAogICAgICAgICAgICBwcmlvciA9IGMoCiAgICAgICAgICAgICAgc2V0X3ByaW9yKCJub3JtYWwoMTAsMTApIixjbGFzcz0iSW50ZXJjZXB0IiksCiAgICAgICAgICAgICAgc2V0X3ByaW9yKCJub3JtYWwoMCwyMCkiLGNsYXNzPSJiIiksCiAgICAgICAgICAgICAgc2V0X3ByaW9yKCJjYXVjaHkoMCwxKSIsY2xhc3M9InNpZ21hIikKICAgICAgICAgICAgKSwKICAgICAgICAgICAgaXRlcj01MDAwLAogICAgICAgICAgICB3YXJtdXA9MTUwMCwKICAgICAgICAgICAgZGF0YT10b21hdG8pCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KaW50MiA8LSB1cGRhdGUoaW50MSwgZm9ybXVsYS4gPSBpbnRsZW5nIH4gYWNzKnRydCArICgxfHNwZWNpZXMpICsgKDF8c2hlbGYpKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CmludDMgPC0gdXBkYXRlKGludDIsIGZvcm11bGEuID0gaW50bGVuZyB+IGFjcyp0cnQgKyAoMXxzcGVjaWVzKSkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyfQppbnQ0IDwtIHVwZGF0ZShpbnQzLCBmb3JtdWxhLiA9IGludGxlbmcgfiBhY3MqdHJ0ICsgKDF8c2hlbGYpKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CmludDUgPC0gdXBkYXRlKGludDEsIGZvcm11bGEuID0gaW50bGVuZyB+IGFjcyp0cnQpCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KbG9vKGludDEsaW50MixpbnQzLGludDQsaW50NSkKYGBgCgpTbyB3ZSBzaG91bGQgZ28gd2l0aCBpbnQ0IChubyBzcGVjaWVzIGVmZmVjdCkuICBDYXNlIGNvdWxkIGFsc28gYmUgbWFkZSBmb3IgaW50NSAobm8gcmFuZG9tIGVmZmVjdHMpCmBgYHtyfQpzdW1tYXJ5KGludDQpCnBsb3QoaW50NCwgYXNrPUZBTFNFKQpyYW5lZihpbnQ0KQpgYGAKCmBgYHtyfQppbnRsZW5nLnByZWRpY3Rpb24gPC0gY2JpbmQocHJlZC5kZixmaXR0ZWQoaW50NCxuZXdkYXRhID0gcHJlZC5kZiwgcmVfZm9ybXVsYSA9IE5BKSkKY29sbmFtZXMoaW50bGVuZy5wcmVkaWN0aW9uKVs2OjddIDwtIGMoImxvdzk1IiwiaGlnaDk1IikKaGVhZChpbnRsZW5nLnByZWRpY3Rpb24pCndyaXRlLmNzdihpbnRsZW5nLnByZWRpY3Rpb24sZmlsZT0iaW50bGVuZ19hY2Nlc3Npb24uY3N2IikKYGBgCgpgYGB7cn0KcGwgPC0gZ2dwbG90KGludGxlbmcucHJlZGljdGlvbixhZXMoeD1hY3MseT1Fc3RpbWF0ZSx5bWluPWxvdzk1LHltYXg9aGlnaDk1LCBmaWxsPXRydCkpCnBsIDwtIHBsICsgZ2VvbV9iYXIocG9zaXRpb249ImRvZGdlIixzdGF0PSJpZGVudGl0eSIpCnBsIDwtIHBsICsgZ2VvbV9lcnJvcmJhcihwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0uOSksIHdpZHRoPTAuNSkKcGwgPC0gcGwgKyBmYWNldF93cmFwKCB+IHNwZWNpZXMsIG5jb2w9Miwgc2NhbGVzPSJmcmVlX3giKQpwbCArIGdndGl0bGUoImludGVybm9kZSIpCmdnc2F2ZSgiaW50ZXJub2Rlc19hY2Nlc3Npb24ucGRmIixoZWlnaHQ9MTEsd2lkdGg9OC41KQpgYGAKCiMjIHNwZWNpZXMgbGV2ZWwgYW5hbHlzZXMKCiMjIyBoeXBvY290eWwKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwNiA8LSBicm0oaHlwIH4gc3BlY2llcyp0cnQgKyAodHJ0fGFjcykgKyAoMXxzaGVsZiksCiAgICAgICAgICAgIHByaW9yID0gYygKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgyMCwxMCkiLGNsYXNzPSJJbnRlcmNlcHQiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoIm5vcm1hbCgwLDEwKSIsY2xhc3M9ImIiKSwKICAgICAgICAgICAgICBzZXRfcHJpb3IoImNhdWNoeSgwLDEpIixjbGFzcz0ic2lnbWEiKQogICAgICAgICAgICApLAogICAgICAgICAgICBpdGVyPTUwMDAsCiAgICAgICAgICAgIHdhcm11cD0xNTAwLAogICAgICAgICAgICBkYXRhPXRvbWF0bykKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwNyA8LSB1cGRhdGUoaHlwNiwgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQgKyAoMXxhY3MpICsgKDF8c2hlbGYpICkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCgpgYGB7ciwgcmVzdWx0cz0naGlkZSd9Cmh5cDggPC0gdXBkYXRlKGh5cDcsIGZvcm11bGEuID0gLiB+IHNwZWNpZXMqdHJ0ICsgKDF8c2hlbGYpICkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwOSA8LSB1cGRhdGUoaHlwNywgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQgKyAoMXxhY3MpICkKc2F2ZS5pbWFnZShmaWxlPSJ+L0dvb2dsZSBEcml2ZS9SZGF0YS9jYXJsb3MuUmRhdGEiKQpgYGAKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaHlwMTAgPC0gdXBkYXRlKGh5cDcsIGZvcm11bGEuID0gLiB+IHNwZWNpZXMqdHJ0KQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3J9CnN1bW1hcnkoaHlwNikKYGBgCgpgYGB7cn0Kc3VtbWFyeShoeXA3KQpgYGAKCmBgYHtyfQpzdW1tYXJ5KGh5cDgpCmBgYAoKYGBge3J9CnN1bW1hcnkoaHlwOSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShoeXAxMCkKYGBgCgpgYGB7cn0KbG9vKGh5cDYsaHlwNyxoeXA4LGh5cDksaHlwMTApCmBgYAoKaHlwNiwgNyBlcXVpdmFsZW50IChubyBuZWVkIGZvciB0cnR8YWNzKQpoeXA4LCAxMCBtdWNoIHdvcnNlIChrZWVwIGFjcyBpbnRlcmNlcHQgcmFuZG9tIGVmZmVjdCkKaHlwOSBhIGJpdCB3b3JzZSAoc29tZSBldmlkZW5jIGZvciBzaGVsZiBlZmZlY3RzKQoKZ28gd2l0aCBoeXA3CgpgYGB7cn0KcHJlZC5kZi5zcGVjaWVzIDwtIHVuaXF1ZSh0b21hdG9bLGMoInNwZWNpZXMiLCJ0cnQiKV0pCmh5cC5zcGVjaWVzLnByZWRpY3Rpb24gPC0gY2JpbmQocHJlZC5kZi5zcGVjaWVzLGZpdHRlZChoeXA3LG5ld2RhdGEgPSBwcmVkLmRmLnNwZWNpZXMsIHJlX2Zvcm11bGEgPSBOQSkpCmNvbG5hbWVzKGh5cC5zcGVjaWVzLnByZWRpY3Rpb24pWzU6Nl0gPC0gYygibG93OTUiLCJoaWdoOTUiKQpoZWFkKGh5cC5zcGVjaWVzLnByZWRpY3Rpb24pCndyaXRlLmNzdihoeXAuc3BlY2llcy5wcmVkaWN0aW9uLGZpbGU9Imh5cF9zcGVjaWVzLmNzdiIpCmBgYAoKYGBge3J9CnBsIDwtIGdncGxvdChoeXAuc3BlY2llcy5wcmVkaWN0aW9uLGFlcyh4PXNwZWNpZXMseT1Fc3RpbWF0ZSx5bWluPWxvdzk1LHltYXg9aGlnaDk1LCBmaWxsPXRydCkpCnBsIDwtIHBsICsgZ2VvbV9iYXIocG9zaXRpb249ImRvZGdlIixzdGF0PSJpZGVudGl0eSIpCnBsIDwtIHBsICsgZ2VvbV9lcnJvcmJhcihwb3NpdGlvbj1wb3NpdGlvbl9kb2RnZSh3aWR0aD0uOSksIHdpZHRoPTAuNSkKcGwgKyBnZ3RpdGxlKCJoeXBvY290eWxzIHBlciBzcGVjaWVzIikKZ2dzYXZlKCJoeXBvY290eWxzX3NwZWNpZXMucGRmIikKYGBgCgojIyMgaW50ZXJub2RlCgpgYGB7ciwgcmVzdWx0cz0naGlkZSd9CmludDYgPC0gYnJtKGludGxlbmcgfiBzcGVjaWVzKnRydCArICh0cnR8YWNzKSArICgxfHNoZWxmKSwKICAgICAgICAgICAgcHJpb3IgPSBjKAogICAgICAgICAgICAgIHNldF9wcmlvcigibm9ybWFsKDIwLDEwKSIsY2xhc3M9IkludGVyY2VwdCIpLAogICAgICAgICAgICAgIHNldF9wcmlvcigibm9ybWFsKDAsMTApIixjbGFzcz0iYiIpLAogICAgICAgICAgICAgIHNldF9wcmlvcigiY2F1Y2h5KDAsMSkiLGNsYXNzPSJzaWdtYSIpCiAgICAgICAgICAgICksCiAgICAgICAgICAgIGl0ZXI9NTAwMCwKICAgICAgICAgICAgd2FybXVwPTE1MDAsCiAgICAgICAgICAgIGRhdGE9dG9tYXRvKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3IsIHJlc3VsdHM9J2hpZGUnfQppbnQ3IDwtIHVwZGF0ZShpbnQ2LCBmb3JtdWxhLiA9IC4gfiBzcGVjaWVzKnRydCArICgxfGFjcykgKyAoMXxzaGVsZikgKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKCmBgYHtyLCByZXN1bHRzPSdoaWRlJ30KaW50OCA8LSB1cGRhdGUoaW50NywgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQgKyAoMXxzaGVsZikgKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3IsIHJlc3VsdHM9J2hpZGUnfQppbnQ5IDwtIHVwZGF0ZShpbnQ3LCBmb3JtdWxhLiA9IC4gfiBzcGVjaWVzKnRydCArICgxfGFjcykgKQpzYXZlLmltYWdlKGZpbGU9In4vR29vZ2xlIERyaXZlL1JkYXRhL2Nhcmxvcy5SZGF0YSIpCmBgYAoKYGBge3IsIHJlc3VsdHM9J2hpZGUnfQppbnQxMCA8LSB1cGRhdGUoaW50NywgZm9ybXVsYS4gPSAuIH4gc3BlY2llcyp0cnQpCnNhdmUuaW1hZ2UoZmlsZT0ifi9Hb29nbGUgRHJpdmUvUmRhdGEvY2FybG9zLlJkYXRhIikKYGBgCgpgYGB7cn0KbG9vKGludDYsaW50NyxpbnQ4LGludDksaW50MTApCmBgYAoKaW50NiBpcyBiZXN0CgpgYGB7cn0KaW50LnNwZWNpZXMucHJlZGljdGlvbiA8LSBjYmluZChwcmVkLmRmLnNwZWNpZXMsZml0dGVkKGludDYsbmV3ZGF0YSA9IHByZWQuZGYuc3BlY2llcywgcmVfZm9ybXVsYSA9IE5BKSkKY29sbmFtZXMoaW50LnNwZWNpZXMucHJlZGljdGlvbilbNTo2XSA8LSBjKCJsb3c5NSIsImhpZ2g5NSIpCmhlYWQoaW50LnNwZWNpZXMucHJlZGljdGlvbikKd3JpdGUuY3N2KGludC5zcGVjaWVzLnByZWRpY3Rpb24sZmlsZT0iaW50X3NwZWNpZXMuY3N2IikKYGBgCgpgYGB7cn0KcGwgPC0gZ2dwbG90KGludC5zcGVjaWVzLnByZWRpY3Rpb24sYWVzKHg9c3BlY2llcyx5PUVzdGltYXRlLHltaW49bG93OTUseW1heD1oaWdoOTUsIGZpbGw9dHJ0KSkKcGwgPC0gcGwgKyBnZW9tX2Jhcihwb3NpdGlvbj0iZG9kZ2UiLHN0YXQ9ImlkZW50aXR5IikKcGwgPC0gcGwgKyBnZW9tX2Vycm9yYmFyKHBvc2l0aW9uPXBvc2l0aW9uX2RvZGdlKHdpZHRoPS45KSwgd2lkdGg9MC41KQpwbCArIGdndGl0bGUoImludGVybm9kZSBwZXIgc3BlY2llcyIpCmdnc2F2ZSgiaW50ZXJub2Rlc19zcGVjaWVzLnBkZiIpCmBgYA==